Enhancement of radiation trapping for quasi-resonant scatterers at low temperature 
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We present a transport equation for the incoherent propagation of radiation inside a quasi-resonant 
atomic gas at low temperature. The derivation is based on a generalized Bethe-Salpeter equation 
taking into account the motion of the atoms. The obtained equation is similar to the radiative 
transfer equation. It is solved numerically by an original Monte Carlo approach in the case of a slab 
geometry. The partial frequency redistribution caused by the small velocity of the scatterers make 
the emitted flux outside the system and the energy density inside the medium to behave differently 
than in the case of complete frequency redistribution. In particular, the long time dependence of 
the specific intensity (escape factor) is slightly different from the Holstein prediction. 

PACS numbers: 42.25.Dd, 32.80.-t 
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I. INTRODUCTION 

The study of light propagation in a scattering system 
has become a very active field of research in mesoscopic 
physics for few years because of its many apphcations in 
particular in biomedical imaging [l| . Lots of key features 
like Anderson localization of photons 0] or fluorescence 
lifetime of a single emitter embedded in a complex sys- 
tem have to be understood in such a system [3,]. A very 
good medium to deal with these fundamental phenomena 
is composed of cold atoms. Cold atoms can be manipu- 
lated with a high degree of control in particular to realize 
a monodisperse ensemble of strongly resonant point scat- 
terers, free of defects and absorption. 

Since the thirties and the pioneering work of Kenty 4] , 
the trapping of photons in a scattering and resonant sys- 
tem is well understood. For a dilute system illuminated 
by a laser beam and in the limit of low saturation of the 
atomic transition (e.g. at low light intensity for resonant 
light) and sufficiently low density for atom-atom coUi- 
sional broadening to be neglected, the scattering remains 
fully elastic. Thus, during its propagation through the 
system, the light remains monochromatic and the scat- 
tering mean-free path can be continuously adjusted by 
tuning the laser frequency. It is then possible to observe 
coherent effect such as weak and strong localization 

At very low temperature and if the system is not too 
large, the residual motion of the atoms does not modify 
significantly the frequency of a scattered photon and co- 
herent effects such as coherent backscattering 0, @j can 
still be observed, as if the atoms were pinned at fixed 
positions. Residual motion even provides us with a sig- 
nificant advantage: averaging over different disorder re- 
alizations is done automatically by just waiting for the 
atoms to move by say one wavelength, i.e. by trivial time 



averaging. 

The Doppler shift of the scattered photon can be ne- 
glected if the frequency change does not modify the prop- 
agation properties, meaning that the scattered photon 
can be rescattered like the incoming ones. The scattering 
cross-section varies rapidly around an atomic resonance 
at frequency wq over a frequency range F, where F is the 
inverse of the lifetime of the atomic excited state. Thus, 
if k denotes the wave factor of the incoming photons and 
V the typical atomic velocity, the typicla Doppler shift 
after a scatteting event is kv and the regime of "quasi- 
pinned" disorder is reached when: 



(1) 



As will be seen in the following, this is only a very rough 
criterion and a more careful analysis is needed. Eq. ^ 
can be rewritten as i;/c <gC wq/F, the latter quantity be- 
ing the quality factor of the atomic resonance, a number 
typically of the order of 10*. This thus puts a severe limit 
on the atomic velocity, of the order of 1 m/s, and explains 
why observation of coherent effects requires cold atomic 
gases. 

However, even if the temperature is low, for a suf- 
ficiently large medium, many scattering events can be 
chained before a photon escapes. Even if the Doppler 
shift of a single scattering event is small, the global 
frequency shift accumulated along a multiple scattering 
path can be non negligible. This depends of course on 
the optical thickness b: 



(2) 
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where L is the size of the medium and I the scattering 
mean free path of the photon in the medium. At large b, 
a typical multiple scattering path can be viewed as a ran- 
dom walk of the photon in the medium with step Thus 
the photon is multiply scattered about b^ times before 
escaping. As the Doppler shift depends on the relative 
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orientation of the atomic velocity and the incoming and 
outgoing wave vectors of the photon, it is on average a 
random variable with zero average and about {kv)"^ vari- 
ance. If successive scattering events are statistically inde- 
pendent, the photon frequency performs itself a random 
walk with step about kv, and the typical accumulated 
Doppler shift after N scattering events is of the order of 
y/Nkv. As TV ~ 6^, a rough criterion for neglecting this 
effect becomes: 

bkv < r. (3) 

Some theoretical works has been done previously in 
particular in rather dense atomic vapors typical of in- 
terstellar atmospheres, when collisional broadening is 
larger than the natural linewidth of the atomic reso- 
nance. In such a regime, the scattering process is in- 
coherent and inelastic: the frequency of the scattered 
photon is essentially decorrelated from the incoming fre- 
quency. This is the complete frequency redistribution 
(CFR) regime 0, H, Q ■ Depending on parameters such as 
the optical thickness and the Doppler broadening, several 
regimes can be obtained, leading to various decay rate 
equations such as the Holstein's equation [l^l- In that 
case, the photon frequency can be Doppler shifted very 
far from resonance, so that the medium becomes almost 
transparent: this is a regime where photons trajectories 
can be Levy flights [rTJ. 

The present work deals with the case of low tempera- 
ture (or slow atoms) for which we cannot assume a com- 
plete frequency redistribution. The partial frequency re- 
distribution (PFR) is then due to Doppler effect. This is 
an interesting new regime where the radiation transport 
has both coherent (the individual scattering event) and 
incoherent (because of randomization due to the atomic 
motion) aspects. Moreover, the study of such a system is 
of primary importance to design new high capacity quan- 
tum memories using cold atoms, in particular to give a 
reliable expression of the escape factor (i.e. exponential 
time decay rate at long times) . The main idea is to derive 
a transport equation for the incoherent radiation from 
first principles in the case of moving atoms, taking into 
account the resonant character of the scatterers. This is 
very different from e.g. diffusing-wave spectroscopy the- 
ory m, or electron transport . It is also expected 
that the dynamics of the atoms affect the coherent prop- 
erties of the system as shown experimentally [l^ . Theo- 
ritical and numerical works are in progress in particular 
to describe the thermal decoherence of the backscattering 
cone. 

More precisely, we deal with systems composed of a 
cloud of two-level atoms: this is a very good approxima- 
tion for quasi-resonant atomic scatterers as e.g. Rubid- 
ium atoms used in experiments. As we consider here only 
the incoherent transport of radiation, the existence of an 
internal hyperfine structure which leads to considerable 
modifications of the coherent transport such as coherent 
back scattering [l^, is irrelevant. We assume a dilute gas 
where the scattering mean free path is much larger than 



the optical wavelength. We will also treat atoms as clas- 
sical point scatterers. This assumption certainly fails at 
very low temperature where the atomic de Broglie wave- 
length becomes comparable to the optical wavelength: 
in this regime, the recoil effect induced by the scattering 
of a single photon significantly affects the motion of the 
center of mass of the atom. We thus exclude ultra-cold 
atomic gases. Properly taking into account the quantum 
nature of the external atomic motion is a much more 
complicated problem, see e.g. [l3| for the simple case of 
two atoms. For usual atoms like Rubidium, the recoil 
velocity is of the order of few mm/s, and one can easily 
have cold atoms faster than the recoil velocity, but still 
obeying Eq. ([T]). We also exclude collective quantum ef- 
fects such as Bose-Einstein condensation, taking place 
when the atomic de Broglie wavelength is comparable to 
the inter-atomic distance. 

We are interested in typically an atomic cloud in a 
magneto-optical trap (MOT), scattering photons from a 
laser beam. We will thus assume a thermal distribution of 
the atomic velocities. In such a system, the radiation felt 
deep in the medium is quasi-isotropic such that atoms 
are not significantly pushed by incoming photons. On 
the other handm, in a real experiment, if the laser beam 
is sufficiently intense, the accumulated recoil effect may 
significantly disturb the atomic velocity distribution. For 
simplicity, we neglect this effect, although taking it into 
account would be easy in the present framework. 

The paper is organized as follows: in Sect. |TT1 we 
start from first principles to obtain a generalized form of 
the well-known Dyson (for the average field) and Bethe- 
Salpeter (for the field correlations, including the aver- 
age intensity) equations. Then, some assumptions are 
done in Sect. ITTTl to derive the transport equation which 
is the main result of this paper. This equation is solved 
numerically in Sect. IIVI using an original Monte Carlo 
scheme. Finally, Sect. |V] is devoted to the modal study 
of the transport equation to explicitly derive the spectral 
and temporal behaviors at large scales. In particular, we 
show that, although the spatial motion of the photons is 
in general not diffusive, it is possible to derive equations 
of the Fokker-Planck type, which govern the evolution of 
the intensity distribution as a function of position and 
frequency. This section ends up with the simple picture 
story of photons migration valid at large scales. 



II. DYSON AND BETHE-SALPETER 
EQUATIONS FOR MOVING SCATTERERS 

A. Scattering operator 

Starting with first principles means that we have to de- 
rive a generalized form of the scattering operator taking 
into account the Doppler effect. It will be the funda- 
mental quantity of our derivation. As mentioned above, 
the recoil is not taken into account in this derivation 
because the recoil-induced drift of the frequency distri- 
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bution is usually much smaller than the Doppler-induced 
spreading (urecoii/'*^ ~ 10"^ — 10~^). Note that going be- 
yond this approximation is far from trivial. We denote by 
t (k, k', (jj) the scattering operator for a single fixed atom, 
k and k' are the incident and scattered wave-vectors re- 
spectively and Lu is the frequency. The reader will find in 
Appendixl^the conventions used for the spatio-temporal 
Fourier transforms of all quantities. For the sake of sim- 
plicity, we assume that all the scattering atoms are iden- 
tical. Let us consider an atom moving at velocity v con- 
stant during the whole scattering process. We denote by 
t^j (k, k', a;, cj') the scattering operator for this moving 
atom (subscript m means moving). Thus, the incident 
(iJinc) and scattered (£^sca) fields are related by 

£^sca (k,C^)= J Go (k,C^)C (k,k',C^,C^') 

^ .d^k'dc^' 

in the scalar approximation (polarization effects ne- 
glected), where Gq is the Green function in the vacuum. 
Note that taking into account polarization effects is a 
straightforward extension of the present analysis. Indeed, 
Doppler effect acts similarly on all polarization compo- 
nents. It is thus sufficient to add the angular dependence 
of the scattering cross-section in all formulas. We drop 
this dependence for simplicity and focus on the main fea- 
tures related to the motion of the scatterers. 

Considering now this expression in the atom frame, we 
have 



Ssca. (k, uj) 



Go (k, w -I- k • v) 



(4) 



X t^^ (k,k',w-|-k- v,w' 
d^k' duj' 



k' 



X ^inc (k',t^') 
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where £ is the field in the atom frame in which the atom 
is fixed. The relation between the incident and the scat- 
tered fields can be written in the form 



f sea (k, Uj) 



I Go (k,C^)i(k,k',C^)finc (k',C 



d^k' 



87r3 ■ 
(5) 

By identifying equations ([4]) and ^ and taking into ac- 
count that the velocity of the atom is weak enough to 
have k • V <C w in the Green operator, we have 



(k, k , L 



--2TTt (k,k',w - 
X 8 \uj' — UJ 



k- v) 

fk' - k) 



(6) 



which is the expression of the generalized scattering oper- 
ator, a fundamental quantity in our derivation. 5 denotes 
the Dirac-delta function. The physical interpretation of 
tm is straightforward: the modification of frequency in 
the t operator corresponds to the Doppler shifted fre- 
quency in the atomic frame, while the 5 function ex- 
presses that the scattering is elastic in the atomic frame. 



The latter expresses a strong correlation between the in- 
coming and the outgoing frequency, in strong contrast 
with the usual complete frequency redistribution hypoth- 
esis 0,i|. 



B. Two-levels atom polarizability 

The scalar scattering operator for pinned point scat- 
terers is related to the polarizability by 

^2 

t (k, k', Lj) = ja {uj) 

where cq is the light velocity in vacuum. The system is 
composed of two- level atoms with resonant frequency Wq. 
r is the spontaneous emission decay rate. The incoming 
photons are quasi-resonant, such that 5 = uj — ujo ujo- 
The two-level atom scalar polarizability is then given by 



a {uj) 



Att r/2 



iT/2 



(7) 



where ko — ujo/cq is the wave- vector in vacuum at fre- 
quency ujo- We thus have 



t(k,k',c^) 



An r/2 
'k^5 + ir/2' 



(8) 



Thus, the generalized scattering operator is obtained by 
plugging Eq. ^ in Eq. ^ . The main differences between 
this equation and the expression of the scattering opera- 
tor for fixed atoms are the Dirac S function representing 
the Doppler frequency shift. 



C. Dyson equation 

The Dyson equation is a closed equation for the field 
(E) averaged over an ensemble of realization of the scat- 
tering medium [l^, [l^. For the case of fixed atoms, the 
average is performed over the positions of the scatterers. 
For moving atoms, the average is performed over the ini- 
tial positions, times and velocities of the scatterers. In 
this case, the generalized form of the Dyson equation is 
given by 

(E) (r, t) - Ei,-,, (r, i) + /gq (r - r', t - t') M„ (r', r", t', t") 



X {E) (r",t")d3r'd3r"di'di", 

(9) 

where Mm is called the mass operator. It essentially 
describes the extinction phenomenon due to scattering. 
In particular, it contains the expression of the scatter- 
ing mean-free path depicting the exponential decay of 
the averaged (or coherent) field. The mass operator 
is the equivalent of the self-energy for the Schrodinger 
equation for matter waves in a random potential [14| . 
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For fixed scatterers, one usually considers a statistically 
translationally invariant infinite medium. In our case, we 
also assume that the atom velocity distribution is time- 
independent and uncorrelated with the position, so that 

(r', r", t', t") = M„, (r' - r", t' - i") . 

This expression allows us to rewrite the Dyson equation 
in term of a Fourier transform: 

(E) (k, iu) = Ein, (k, + Go (k, Lu) Mra (k, iu) (E) (k, . 

(10) 

Writing the same equation for the averaged Green func- 
tion, we obtain 

(G) (k, io) - Go (k, ij) + Go (k, ij) M,„ (k, Lo) (G) (k, lu) 



which gives the expression 
(G) ik,u;) = 



2/c2_p_M„(k,c^) 



(11) 



considering that Go (k, w) = 1/ (uj'^/cq — k^). Then we 
can define a wave-vector in the effective homogeneous 
medium k^g = w^/cq — Mm (koff,w), the real and imag- 
inary part of which give the effective wavelength and 
the scattering mean-free path respectively. To derive 
the expression of the mass operator we assume that 
the system is dilute enough for the first order diagram- 
matic expansion of the operator to be valid. This is the 
Foldy-Twersky approximation (or the Born approxima- 
tion) [1^. In that case, the mass operator is the sum over 
all the N scatterers of the averaged scattering operator 
for all accessible positions, times and velocities: 

N 



i=l 



where P (r, t, v) is the probability density to have a 
scatterer of velocity v at position r and at time t. 
For our translationally invariant medium, P (r, t, v) = 
g (v) / (VT) where V is the volume of the system and 
T the time window. In the large V, T limit, the mass 
operator in the Fourier domain reads 



M™ (k, u) 



lim ^ / C (k, k, w, w) 5 (v) dv 

r^oo 1 J 



where p = N/V is the density of scatterers. This leads 
to the final expression of the mass operator if we remark 
that limT^oo Stt^ (lu — lu) /T — 1: 



Mm {k,u;) = p J t(k,k,u;-k-v)g (v) dv. 



(13) 



(12) 



For a thermal distribution, this gives a standard Voigt 
profile for absorption [l^l , which itself reduces to a usual 
Doppler Gaussian profile in the limit kv ^ F. In the lat- 
ter case, the averaging over the atomic velocity makes 
the correlation between the incoming and scattered fre- 
quencies much smaller than for slow atoms. 

D. Bethe-salpeter equation 

The Bethe-Salpeter equation is a closed equation for 
the field autocorrelation function (EE*) [H^, ill. It 
contains an operator K depending on four space variables 
called the intensity (or vertex) operator and describing 
the correlation between two scattering processes. As for 
the Dyson equation, the idea is to generalize this equation 
to the case of moving scatterers. This leads to a new 
vertex operator Km depending on four space variables 
and four time variables. We assume that no source is 
present (the incident current densities associated with 
the coherent term are missing). Thus the generalized 
Bethe-Salpeter equation writes: 



(r'l , t ; ) i;* (r^ , ) - y (G (r'l , r 1 , t'l , 1 1 ) ) (G* (r^ , r2 , , i2 ) ) if™ (r 1 , ra , r2 , r4 , 1 1 , t3 , t2 , t4 ) 

X {E (ra, is) E* {u, U)) A^Yid?Y2d?Y^A^Yi Ahdt2dUAU, (14) 

where the symbol * denotes the conjugate quantity. The global translational invariance in space and time for the 
averaged quantities leads to a factorization of the vertex operator as follows 

Km (ki,k3,k2,k4, wi, W3, W2, a;4) =87r^<5 (ki - ka - ka + k4) x 2ti5 [uji - - ^2 + W4) 

X Km (ki,k3,k2,k4,wi,W3,W2,a;4) . (15) 
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This leads to the following expression of the Fourier transform of the Bethe-Salpeter equation for an infinite medium: 



X Kr, 



,,q,/,q, q,/ q , ^ f , ^ ^ , ^ 

kH — . k H — ,k ,k , cij H ,uj H .uj 

2' 2 2 2' 2' 2' 2 2 



/(k',q,c.',f]) 



d=^k' dw' 



87r3 27r 



(16) 



where the function / is the Fourier form of the spatio-temporal correlation function of the electric field (i.e. 
/ (k, q, , c^, f]) = (i; (k + q/2, + f^/2) i;* (k - q/2, - 17/2))). 

As for the mass operator, we only keep the first term of the diagrammatic expansion of the vertex operator (Ladder 
approximation valid for a dilute gas) which leads to a sum over all scatterers, and an averaging over all accessible 
positions, times and velocities of the scattering operator correlation function. In Fourier space, we obtain: 



ifm (ki,k3,k2,k4,cji,Lj3,cj2,t^4) = lim — / C (ki, ka, t^i, W3) (kj, k4, c^a, t^4) 5 (v) dv. 



The averaging over the atomic velocity is performed over the product t^ji^* to take into account the fact that the 
scattering process occurs on the same atom for the field and its conjugate. Using this expression of the scattering 
operator in Eq. ([5]), the intensity operator becomes 



,,q,/,q, q,/ q , ^ i , ^ ^ , ^ \ i . 

kH — ,k H — ,k ,k ,cijH ,0; H ,0; ,uj p / t 

2 2 2 2 2 2 2 2 ' ' 



k+77:k +-,tj+--(k+-)-v 



X t* 



1 ^ 1 ' ^ 

k , k ,uj - 

2' 2' 



2 



2' 2' 2 

27:6 [lj' -Lu- (k' - k) • v] g (v) dv 



(17) 



III. TRANSPORT EQUATION 



The transport equation we will obtain is an equation governing the specific intensity / (r, u, t, lo) inside the system. 
This is a local (r) and directional (u) radiative flux at time t and frequency uj. This quantity will be defined using 
the field autocorrelation function. So the root of the derivation is the Bethe-Salpeter equation in the form obtained 
in Eq. (|16p . To exhibit the spatio-temporal derivatives of the specific intensity, we can transform the product of the 
averaged Green functions in a difference: 1/{AB) — {l/A — 1/ B)/{B — A), and use the averaged Green function in 
Eq. ini)- We obtain: 

(K^-!-§))(-(-f-§)> 



G* 



-, UJ - 



— — + 2k • q + M„ 



, q ^ 

'^+2'"+2 



M* 



q n 

k — — , — — 
2' 2 



(18) 



The effective wave-vector in the medium is given by 
aj^/c§ — Mm (kcff, w). By replacing the mass op- 
erator by its approximate expression, Eq. (|13p and the 
polarizability by Eq. ([7]), we obtain 



3 



cl ka 



r/2 ((5-k-v) 



{S-k-v) 

r2/4 



(J - k • v)-" + rV4 



-r2/4 

g (v) dv. 



g (v) dv. 



In a dilute medium such that pXq ^ 1 where Aq = 27r/fco 
is the wavelength in the vacuum at frequency ujq, we have 
3 < 3? ^ uj^cl Writing k,s = k' + ik" we 



find 



^cffj 



27r a; 

Aoff Co ■ 

1 



(19) 



2^ 



where Acff is the effective wavelength (close to the vac- 
uum wavelength in a dilute medium) and £ ^ Xq the 
scattering mean-free path. Physically, this means that 
the spatial variations of the specific intensity take place 
on a scale much longer than the wavelength. We can 
thus restrict ourselves to g <C fc, k' in Eq. (|16p . Similarly, 
the temporal variations of the specific intensity are slow 
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(typically of the order of T~^) compared to the tempo- 
ral variations of the wave (on a time scale oJq^). This 
allows to define properly the specific intensity with two 
independent spatial variables (position r and direction 
u) and two independent temporal variables (time t and 
frequency lu) and assume <ti lu in Eq. ()16p . Note how- 



ever that we do NOT assume <C F — meaning that we 
keep the full dynamics on the time scale — and that 
because the polarizability varies around the resonance on 
a scale F, we must keep explicitly this dependence. Us- 
ing these approximations coupled to Eqs. p^ . (fTT)) and 
(fT8|) . the Bethe-Salpeter equation (fTH]) reduces to 



2nu. 



2k<i-p a k,k,c^-k.v + - k,k,c^-k.v-- 5(v)dv 



n 



X p J t ^k, k', — k • V 



i* k,k',w-k- V 



87r3 ■ 
(20) 



(k' - k) • v] g (v) / (k', q, uj' , ft) dw'dv 



We can glimpse in this equation the role of each term in 
the coming transport equation. The first two ones corre- 
spond respectively to the temporal and spatial evolution 
of the specific intensity. The term t — t* is related to the 
extinction coefficient and the integral on the right side 
of the equation with the term tt* refers to the scattering 
process. 

In a dilute system, the imaginary part of the av- 
eraged Green function is peaked and can be written 
(this is the "on-shell" approximation) Q[(G(k, w))] w 



— 7r<5 [w^/cq — k^] . As w is closed to wq, we can replace 
k and k' by /cq- Finally, defining the specific intensity / 
as the Fourier transform of the autocorrelation function 
of the field (i.e. the so-called Wigner transform of the 
field) [H 

/ (k, q, LO,n)^S (fco -k)I (u, q, S, n) (21) 

where we recall that S = lu — ujq, we obtain the Fourier 
form of the transport equation 



in 

h m • q 

Co 



fie (S - fcou ■v,n)g (v) dv 



I{u,q,S,n)^^ [ f [ ps{S-koU-^,n)d[S' -S-ko{u' -u)■^r] 

X g (v) / (u', q, S', Q) dvdu'd(5' (22) 



where 



fie [S - koU ■ V, fl) 



ip_ 
2ko 
P 



t { LU ~ koU ■ V 



and fis {S — koU • v, f2) = —t ( lu — koU • v + — ) t* ( a; — /cqU • v 



t* [ LU — kgU ■ V 



n 



(23) 



where we have written t (k, k', w) = t (lu) for the sake of simplicity. The coefficients fie and fis are like extinction and 
scattering coefficients respectively but depending on time and frequency. The form of the previous equation in real 
space is 



Co dt 



U • Vr 



poo 

I (u, r,5,t) = - / fie {t', d-koU-v)g (v) / (u, r,5,t- t') dvdt' 

Jt'=0 

+ _/ / / / ii{t',S~koU-v)S[5' -S-ko{u' -u}-v]g{v)I{u',r,t-t',S')dvdu'dS'dt' (24) 



where fie and fis are the inverse temporal Fourier trans- 
forms of fie and fis respectively. 



In the following, we suppose that we have a Maxwell- 
Boltzmann distribution of velocities given by 



3(v) 



1 



[vV2^' 



■ exp 



V 

2^ 



(25) 



7 



that is a gaussian distribution with vanishing mean and 
standard deviation v. Eq. ([M)) is the main result of this 
paper. It is vaUd whatever v. But in the following, we 
focus on kov ^ F (Eq. It has the same structure 

as the well-known radiative transfer equation (RTE) as 
derived by Chandrasekhar [2^ except that it exhibits a 
frequency coupling and temporal convolution products. 
Note also that the phase function p (u, u') (part of an 
incident beam in the direction u' scattered into the di- 
rection u) is constant and equal to l/iir because of the 
isotropic scattering assumption. As already mentioned, 
polarization effects can be easily taken into account by 
inserting in Eq. p4)) the full phase function of e.g. a 
Rayleigh or a Mie scatterer. 

The coupling between the components at various de- 
tunings S is due to the motion of the scatterers: the 
Doppler shift at each scattering event produces a change 
of frequency for the scattered photon. For a pinned sys- 
tem, V ^ 0, the Gaussian becomes a Dirac-delta func- 
tion, and the components at various detunings S become 
uncoupled. For each component, one recovers the well- 
known RTE. 

The convolution products in time describes the fact 
that resonances can be responsible for a drastic reduc- 
tion of the velocity of energy propagation in the medium. 
This is developed in Sect. |Vl Far from resonance, the 
extinction and scattering coefficients do not depend on 
frequency and, again, we recover the well-known RTE. 

The optical theorem (or the Ward identity) is verified 
in Eq. ([M)) . Indeed, we have 



L 



/ie (^ - fcou • V, n) =-5- / / / ^s{S - kau-v,n) 

X5 (v) dv X S[5' — S — ka (u' — u) • v] 

X g (v) dvdu'd^'. 



No absorption is present in the system: the extinction 
phenomena is only due to scattering. More precisely, the 
photon flux is conserved but not exactly the energy flux. 
This is because we neglected the recoil effect: when the 
frequency of the scattered photon is different from the 
incoming frequency, the energy difference is in fact trans- 
ferred to the atomic scatterer whose velocity is slightly 
modified. In our case, the violation of energy conserva- 
tion is of the order of F/ojq and can be safely forgotten. 



The steady-state transport equation for pinned atoms 
makes it possible to define the scattering mean-free path 



£ (S) = lim = ,^ ^, 
n^O fie (S, ft) 



T2" 



(26) 



Losses 



Gains 



where £q = kg / (47rp) is the scattering mean- free path at 
resonance. 

IV. NUMERICAL SIMULATIONS 



A Monte-Carlo scheme 



By integrating over the frequency S' , the structure of 
Eq. is similar to the usual RTE. We have 



in 

h iu • q 

Co 



j fie (S - kou ■ V, n) g (v) dv / (u, q, S, n) =^ J J fi^ {5 - kou ■ v. 



n)g{^) 



X / (u', q, (5 + fco (u' - u) • V, n) dvdu'. 



(27) 



r 



This specific form explains why it is possible to solve it 
using a Monte-Carlo scheme, the only complication be- 
ing that we have to deal with complex probability density. 
The details of the computation are given in appendix IB] 
In a word, our original approach consists in solving the 
temporal Fourier transform of the transport equation and 
not directly the transport equation itself. This allows us 
to avoid the problem of the temporal convolution prod- 
ucts. This Monte-Carlo scheme is the only one which 



we could find to work in this case and should be used 
whenever resonances are present in the system. 



B. Slab geometry 

We numerically study a cold atomic medium whose 
shape is a slab, infinite along the x and y directions, but 
extending along the z direction in the range < z < L. 
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Parameters and numerical values 
i^o/27r = 3.85 x 10'" ^ 



V = 9.17 X 10"^m.s 



Hz 

■ 1 



Z2_ 



1.00 X 10"^ s 



r/27r = 5.9 X 10" Hz 
p = 2.12 X 10^^ m=^ 
L = 40£o or L = lO^o 



TABLE I: Numerical values of the parameters which are cho- 
sen such that all assumptions made to derive the transport 
equation are fulfilled in particular pAp ^ 1, kov <C F and 
T <^ uiQ. These parameters correspond to the Rubidium atom. 



This system is chosen for its simplicity and does not 
represent a real atomic cloud released from a standard 
magneto-optical trap. Nevertheless, this is a true three 
dimensional geometry. The slab is illuminated at nor- 
mal incidence from the left by a gaussian pulse of tem- 
poral FWHM Ti/2 at frequency lol, detuned from the 
atomic resonance hy Sl = ojl — ojq <C ojq. We have 
Ti/2 S> 2-k/u)l so that the incident pulse is well-defined 
and its temporal and spectral evolutions are indepen- 
dent. We are interested in the expression of the energy 
density and the reflected and transmitted fluxes denoted 
by U {z,t,6), R{t,6) and T{t,S) respectively. The sys- 
tem and the notations are clarified in Fig. [T] and all the 
numerical parameters are given in Table [H The compu- 
tations are done in a cluster of 96 Intel Xeon processors 
at 3.0 GHz. 

The most important parameter is the ratio of the 
Doppler shift to the natural linewidth. As discussed 
above, we choose it much smaller than unity: 



kpv 

r 



1 

50 



The optical thickness at resonance is chosen as: 



L 

To 



10 or 40 



(28) 



(29) 



The first value is chosen so that Eq. ^ is satisfied so 
that the effect of the atomic velocity is expected to be 
small and the temporal behavior similar to the one for 
pinned atoms. In the second case, important effects are 
expected. The atomic density is chosen such that: 



pXl = O.K 1 



(30) 



so that the medium is dilute and its index of refraction 
very close to unity. 

This specific choice of parameters for the calculation 
is such that the results depend only on the physically 
important parameters: kov/T, b and S/T. 



Ti/2 

A 



• • • 

• • • rtu) 

L 

FIG. 1: (Color online). Geometry of the system studied nu- 
merically. The transmitted and reflected fluxes are integrated 
over all directions. 
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FIG. 2: (Color online). Spectral distribution of the trans- 
mitted photons for five different times in the case of a thin 
slab of optical thickness b = 10. The spectral broadening is 
due to the accumulation along a multiple scattering path of 
small Doppler shifts (here typically 0.02 F at each scattering 
event). At short time, the broadening is small compared to 
the natural linewidth F, resulting in a Gaussian lineshape. 
At long times, the spectral distribution broadening saturates, 
with an approximate Gaussian lineshape. Note that the total 
transmitted flux is larger at Tt = 20 than at Ft = 10, because 
the photons has to be multiply scattered before escaping in 
the forward direction. The decay at longer time is due to the 
finite trapping time in the medium. 



C. Numerical results 

1. Quadratic regime 

When a system with sufficiently small optical thick- 
ness 6 = 10 is illuminated at resonance (i.e. Sl = 0), we 



expect the photon to escape before the frequency broad- 
ening induced by Doppler effect becomes important. The 
result is shown in Fig.dl At short time, one observes dis- 
tinctly a Gaussian lineshape with linewidth <^ T. This is 
the regime where the photon frequency performs a ran- 
dom walk, and thus a global diffusive behaviour is seen 
in frequency space. Most photons escape the medium 
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FIG. 3: (Color online). Transmission T(t,S) versus time for three different frequencies (left plot) and transmission versus 
frequency for four different times (right plot). We have also plotted the temporal behavior of the transmission for a pinned 
system (i.e. v = 0). The incoming pulse is at resonance 5l = 0, the optical thickness of the medium is fo = 40 and kv/T = 0.02. 
Note that the time scale which governs the temporal evolution is large compared to F"^ while the spectral scale is of the order 
of F~^. This validates the independence of the temporal and spectral variables in the specific intensity. 



before their frequency is significantly shifted. At longer 
time, the frequency dependence of the scattering mean- 
free path cannot be ignored. The width of the spectral 
distribution tends to saturate and the lineshape flattens 
at the center. This is the quadratic regime where the 
decay rate r is quadratic in L as for pinned atoms [l^l . 



2. Doppler regime 

We now study a larger system — optical thickness 
& = 40 — for which the Doppler shift is important. The 
numerical results for the transmitted flux at various fre- 
quencies are given in Fig. [31 The temporal behavior is 
an exponential decay at long times. Importantly, the de- 
cay rate is identical r — 262r^^ for all frequency compo- 
nents, although each component has a different scattering 
mean free path in the medium and consequently a differ- 
ent optical thickness (it varies by a factor 3.5 between 
6 = and 6 = O.ST.) However, the time for this regime 
to settle depends on the frequency. More precisely, since 
many scattering events are needed to create photons in 
the system with large detunings, this time increases with 
6. This regime is studied in details in Sect.|Vl The spec- 
tral behavior exhibits two peaks which is very surprising 
at first glance. This comes from the interplay between 
the variation of the scattering mean-free path with the 
frequency (characterized by T) and the frequency redis- 
tribution due to the Doppler shift (characterized by the 
width of the Gaussian shape v). More precisely, since 
the mean free-path increases with the detuning, for ex- 
ample, £{S = F) = 5^0, a resonant photon is quite ef- 
ficiently trapped, whereas a far detuned photon escapes 
rapidly. Therefore, resonant photons undergo a lot of 



Doppler shifts, populating the other frequencies. On the 
contrary, far detuned photons most likely escape without 
any frequency shift. At a certain frequency, these two ef- 
fects are balancing each other, giving rise to a maximum 
photon population, hence, a two peak structure in the 
transmission. Note that the symmetry observed in the 
spectrum of the transmitted flux is due to the fact that 
the illumination is at 6l = and all the properties of the 
system are even functions of S. This would be no longer 
exactly true if recoil effects were taken into account. 

The time dependence of the energy density at the cen- 
ter of the slab (i.e. for z = L/2) plotted in Fig.|4]is similar 
to the one of the transmitted flux. In particular, the ex- 
ponential decay at long times is the same. Concerning 
the spectral dependence, we also see two peaks which are 
nevertheless attenuated at very long times. The reason 
is the same as for the transmitted flux: trapping of pho- 
tons is more efficient at resonance. The spectrum at long 
times is also studied in detail in Sect. |V] using a modal 
approach. 

If we continue to increase the optical thickness of 
the system, the exponential decay becomes lower (i.e. 
T — > -|-oo when L — > -|-oo) and the spacing between the 
peaks increases. In any case, the decay rate is very dif- 
ferent from the one r (x (L/^o)^r~^ for pinned atoms 
and also different from the one deduced from the Hol- 
stein equation (i.e. r oc (L/£o)\/log [L/{2io)]T-^ Bill)- 
We call this the Doppler regime. 

If the temperature of the cold atomic gas is increased 
(at fixed optical thickness), the temporal decay becomes 
faster (i.e. t — > when v +oo) and the spacing be- 
tween the two peaks in the spectral distribution increases. 

Figure [5] gives an example of calculation for a detuned 
excitation at 6l = F. One observes at short time a 
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FIG. 4: (Color online). Density of photons U {L/2, t, S) at the center of the slab versus time for three different frequencies (left 
plot) and energy density versus frequency for fiour different times (right plot). We have also plot the temporal behavior of the 
transmission for a pinned system (i.e. t; = 0). The optical thickness of the system is 6 = 40. 




FIG. 5: (Color online). Relative transmission T{t,5) versus 
frequency for four different times in the case of a detuned 
excitation Sl = T. The spectral behavior at long time tends to 
be symmetric, as in the resonant case. The optical thickness 
of the system at resonance is 40. 



V. TEMPORAL AND SPECTRAL BEHAVIOR 
AT LONG TIMES 



Analysing analytically the temporal and spectral shape 
at long times and for large systems usually ends up with 
the derivation of a diffusion equation which describes the 
evolution of the energy density. The main advantage of 
a transport equation compared to a dilfusion equation 
is that it contains all transport regimes of photons from 
ballistic to diffusive Despite these drawbacks, the 

diffusion approximation is widely used, in particular in 
biological imaging [2^ [s^, HH because of its simplicity. 
In this section, we first derive the diffusion approxima- 
tion from the transport equation for the cases of pinned 
and moving atoms. To achieve this goal, we use a modal 
decomposition of the specific intensity to obtain a disper- 
sion relation which leads to the expression of the diffusion 
coefficient. 



broadening of the spectral distribution. Soon, this broad- 
ening becomes asymmetric, photons being "attracted" 
towards resonance, simply because there are more effi- 
ciently trapped near resonance. In addition, the height 
of the peak on the other side of the resonance increases. 
Finally, at long times, the spectral shape becomes iden- 
tical to the one of Fig. 21 independently of the initial 
detuning: there, multiple Doppler shifts have erased the 
memory of the initial spectral distribution. 



A. Fokker-Planck equation 



As mentioned above, the residual Doppler broadening 
is usually much smaller (by about two orders of magni- 
tude) than the natural width F in the MOT regime, so 
that Eq. IT]) is valid. Thus, we can perform a second order 
Taylor expansion in kv/T of the extinction and scatter- 
ing coefficients, and of the specific intensity. We cannot 
restrict the expansion to the first order because of the 



symmetry of .g (v). This gives 

fle,s {S - koU ■ V, Q,) ^lle,s ^) - ^qU • v^'^ ((5, 9) 

I (u, q,S — koU ■ V, il) ~/ (u, q, S, D.) — fcou • v/' (u, q, S, il) 
I 



i: 



For the sake of simphcity, we consider the slab geome- 
try mentioned in Sect. IIVI This choice does not reduce 
the generaHty of the derivation, the system being a true 
three dimensional geometry. This permits a simplifica- 
tion of the Fokker-Planck equation given in Eq. ((3T|) by 
integrating over x and y (invariance of the problem by 
translation along x and y) and tp (invariance of the prob- 
lem by rotation around z) where 9 and ip are the usual 
spherical angles. To obtain an analytical derivation, the 
terms u- (u' — u) and ||u — u'|| have to be approximated. 
We denote by cos 9 = u • u' the cosine of the scattering 



where the symbols ' and " denote d/dS and d'^/dS'^ re- 
spectively. By integrating over the velocity of the atoms, 
the transport equation (i.e. Eq. ([27|) ) becomes a Fokker- 
Planck type equation which reads in the Fourier space 



I 

angle. So we have 

u • (u' - u) = cose - 1 and ||u - u'|| = ^2(1 - cos 6). 

For large systems and at long times, many scattering 
events occur. So we can replace the cosine of the scatter- 
ing angle by its average given by the so-called anisotropy 
factor. In our case, the scattering being isotropic (the 
phase function is constant), the anisotropy factor van- 
ishes. Thus we write u-(u' — u) ^ — 1 and ||u — u'|| ^ \f2 
and the Fokker-Planck equation reduces to 



-+zu.q + Me (<5,f7) + ^^Me('5,f^) 



/(u,q,,5,r!) 



47r 



pi, {S, n) I (u', q,5,n)-u- (u' - u) (fco^^)' 



{kovY 



p':{d,n)i{u',(i,d,n) + 



u - u'f {kovY 



du'. (31) 



in 



\-ifiq + fie{S,n) , 

Co 2 



Me [0, n) 



I{fi,q,S,n) = - 



+1 



Ps {S, n) I iix', q, 6, n) + {kovY (<5, ^) I' (m', q, 6, ft) 



(-5, il) I in', q, S, n) + {kovY fis {s, n) I" (//, g, <5, Vt) 



(32) 



where p, = cos 9 is the direction. This equation can be 
simplified assuming that we are looking for large systems 
such that L :$> £q (i.e. g < tt/^q and il < F). At de- 
tuning S, the propagating time between two scattering 
events is of the order of Tpj-op = /cq. On the other 
hand, the Wigner time delay characteristic of the scat- 
tering process is about Tgca = F^^. Then we can assume 
that Tpiop ^ Tsca and consider that the speed of light 



in vacuum is infinite. In other words, the propagation 
of energy in the system is hugely reduced by radiation 
trapping — an effect enhanced by the resonant charac- 
ter of the interaction with the atomic scatterers — so 
that free propagation in the vacuum can be considered 
as instantaneous. Using these approximations, Eq. (|32p 
rewrites 



ifj.q 



/ (/i, q, S,n)^^ [/3 {6) I q, 6,n)+j {S) I' ifi', g, (5, 17) + ry {6) I" (/.', g, S, il)] dfi' (33) 
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where we have defined 



1 



1 



1 



£{5)J ioii + Asyr^) 

2 r 



1-4 



~/ (i + 452/r2)2 



kov 



r y £o(i + 4<52/r2)^ 



ry(<S) 



{kovY 



(kov) 



£(6) 4(i + 4<52/r2)- 

Eq ([55)1 reads in the real space 



9z 



Tdt 



I{ti,z,5,t)^-(3(5)I{ii,z,5,t) 



1 



/? (,5) / z, 5,t) + -i {5) —I iiji', z, <5, t) + 77 {5) g^I (/i', z, S, t) 



d/i'. 



(34) 



The first two terms on the left hand side of Eq. ((34l) rep- 
resent the spatio-temporal evolution of the specific inten- 
sity. The first term on the right hand side describes losses 
due to scattering. The integral term over the direction 
and the second derivative in frequency represent gains by 
scattering in the real and frequency spaces respectively. 
Finally the first derivative in frequency is a consequence 
of a drift in the Doppler process which tends to shift the 
frequency towards the resonance. 



B. Modal approach 

We now expand the specific intensity in the form 

/-|-oo 
^ 9q.s (m, 5) exp [iqz + s{q)t]—. 
-oo -Tit ■ 

(35) 



<q) 



To obtain this expression, we have performed a spatial 
Fourier transform of the specific intensity. For each g, 
there exist several eigenvalues s (q) , which may be com- 
plex or real, whose spectrum may be discrete or contin- 
uous, gq^s is the associated eigenvector which depends 
on direction and frequency only. This decomposition is 
a generalized form of the one used to derive the diffu- 
sion approximation starting from the standard RTE [28| . 
For a system of size L, the dominant q is given by 
q = n/ (L + 2zo) , where zq is the so-called extrapolation 
length [32] , of the order of the mean- free path, account- 
ing for the effects at the interface between the medium 
and the vacuum. For an optically thick medium, 6^1, 
one has zq «C L. Therefore, the mode surviving at large 
time, i.e. having the longest lifetime in the system, cor- 
responds to the eigenvalue sq (t/L) with the lowest (in 
magnitude) negative real part. Inserting Eq. (|35p into 



Eq. dSS]) leads to 



+1 



+7('5)C(A^',.5) + r,(J)<,(/.',<5)]dM'. (36) 
Thus the eigenvector is given by 



9q.s (m, S) 



1 



-1-1 



2[iq^i+il + s/T)P{6)] j_, 
+l{S)g'q^Af^',S)+v{6)g'lAf^',S)] d^'- 

By integrating over the direction, /i, we obtain 

1 



[/3(5)<?,,,(m',<5) 



■ arctan 



{l + s/T)PiS) 



'l{S)9q,s {S) + v{S)9'ls{S)]- 



(37) 



This linear second order differential equation can be seen 
as a generalized Schrodinger problem. Note that this is 
not exactly an eigenvalue problem since the eigenvalue is 
not in front of the associated eigenvector. 



C. Pinned atoms 

For pinned atoms, v is equal to zero. Eq. ([37]) de- 
couples in an independent equation for each frequency 
(5, with the following dispersion relation for the lowest 
eigenvalue sq: 



1 



qi{S) 



arctan 



q£iS) 



i + so/r 



= 1 



In an experiment, the relevant S value is fixed by the 
frequency of the excitation ■ For sufficiently large sys- 
tems — larger than the mean-free path at frequency S, 
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see Eq. ([26)1 — we have q(, [5) ^ 1 which corresponds to 
the diffusion approximation regime. A Taylor expansion 
in the dispersion relation leads to 



So 



fe {5) r 



The , i.e. 1/i^, dependence is a characteristic of the 
quadratic diffusion regime [131 • This leads to the defi- 
nition of the diffusion coefficient V and the velocity of 
energy transport ctr by 



V 



So 

«2 



Ctr {5)i{5) 



with 



Ctr = ^ ((5) r = Co 



L3_ 



1 



4^ 

T2" 



(38) 



The first expression of ctr has a simple physical interpre- 
tation: a single scattering event followed by propagation 
on a mean- free path takes time T~^. Note that this is true 
whatever the detuning is, a highly non-trivial property: 
it turns out that both the Wigner time delay at scattering 
and the time delay during propagation (at the group ve- 
locity) depend on the detuning, but not their sum. This 
property has been used in [33| for simple Monte-Carlo 
simulations. The present, more rigorous approach, justi- 
fies these simple simulations. 

The second expression for ctr is simply the product of 
the vacuum light velocity co by the inverse of the qual- 
ity factor of the atomic resonance, by (up to a numerical 
factor) the number of atoms per cubic wavelength, and 
finally by a factor being unity at resonance. Because the 
quality factor ujq/T of an atomic resonance can be very 
high, typically 10®, the velocity can be strongly dimin- 
ished even in a dilute medium where p/fcg ^ 1. In typical 
experimental situations, it can be reduced by 4 orders of 
magnitude [33]. This is the main characteristic of radi- 
ation trapping by resonant scatterers. These results are 
in perfect agreement with literature [s^, [H, [H, [13, [11] . 

The temporal exponential decay rate at long times is 
then given by 



1 

So 



3L' 



a result also in complete agreement with literature 
[33I ]. In a word, defining a diffusion coefficient for fixed 
scatterers is possible and is even well-understood. The 
atomic nature of the scatterers induces a strong decrease 
of the velocity of energy propagation inside the system. 
This property could be used to design new high capacity 
quantum memories. 



D. Moving atoms 

In the case of moving atoms, the detuning S may range 
from —00 to -l-oo., such that the inequality q£oS'^ /T'^ <C 1 
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FIG. 6: (Color online). Temporal decay rate r of the outgoing 
flux and the energy density versus the system size. Crosses 
correspond to the full Monte-Carlo simulations (i.e. resolu- 
tion of Eq. (|27p ). the solid line is related to the numerical 
resolution of the dispersion relation (Eq. (|37|l 'l. For small L 
(but large enough to be in the diffusive regime), the behavior 
is quadratic as for pinned scatterers (dotted line). This is the 
quadratic regime, which extends roughly up to L/^o ~ 30, as 
predicted by Eq. (|40|l . For very large systems, the behavior 
is completely different from the Holstein prediction (dashed 
line). This is the Doppler regime. The agreement between 
the two numerical methods is very good, the small difference 
coming from approximating the cosine of the scattering angle 
by its average (anisotropy factor). 



cannot hold for all (5, ruling out an approach based on a 
Taylor expansion. Physically, it means that for large de- 
tunings, the mean-free path is larger than the system 
size, invalidating the diffusion approximation. As a con- 
sequence, it is not possible to define a global diffusion 
constant in such a system and the evolution of the en- 
ergy density cannot be described with a diffusion equa- 
tion. The same situation occurs for the CFR case and 
was underlined in 0, [^ . This explains in particular the 
necessity to write an integro-differential transport equa- 
tion [3. 

To the best of our knowledge, it is not possible to ob- 
tain an analytical form of the eigenvalues s from Eq. (j37p . 
In order to make some comparisons with the results given 
by the numerical resolution of Eq. (|27p . we have nu- 
merically solved Eq. ([57)1 by a straightforward shooting 
technique. The results are shown in Fig. [6] The agree- 
ment is excellent, only small deviations being visible for 
extremely large systems. Note that the results deviate 
very significantly from the quadratic behaviour and from 
the Holstein prediction using the CFR assumption. This 
proves that resonant cold atomic media are fundamen- 
tally different from usual hot gases as far as radiation 
trapping is concerned. 
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1. "Diffusion equation" 



2. Energy density and outgoing fluxes 



For large systems such that the optical thickness at 
resonance is large, qio <C 1 and if we look at times longer 
than r~^, we can use directly on Eq. the same ap- 
proximation used for deriving Eq. ([33]) to perform the 
integration over t' . We obtain the following simplified 
transport equation: 



U • Vr 

1 

47r 



ld_\ r g (v) dv 
fdi) J^l{5- kou ■ v) 
.9(v) 



/(u, r, S,t) 



471 Jv - koU ■ v) 
X dvdu' 



I{u',r,S + ko (u'-u) •v,i) 

(39) 



Thus this transport equation is what we can call the "dif- 
fusion equation" for the system at finite temperature. A 
simple physical picture can then be extracted from this 
approximation in term of a Monte-Carlo scheme. Note 
that Eq. is valid whatever v. It is exactly the same 
as the standard time-dependent RTE except that the 
Doppler broadening is present and the energy velocity 
depends on frequency. The Monte-Carlo method is then 
the same. When a photon enters the system at the fre- 
quency Sl, it propagates over a distance whose average 
is 1/ [J^ g (v) d\-/£ {S^ — koU ■ v)] (this is the first term 
on the right hand side). Then it is scattered isotropically 
(the last term with the double integral) and undergoes 
a Doppler drift (from 5 -I- fco (u' — u) • v to 6). Next, it 
propagates again with a new scattering mean-free path. 
Each scattering process takes about T~^. The implemen- 
tation of this method is easier than the one detailed in 
Sect. [Bland gives reliable results at large scales only. 

This simple picture can be used to find a quantita- 
tive criterion to discriminate between the Doppler and 
the quadratic regimes too. As explained in the introduc- 
tion, at large 5, a typical multiple scattering path can be 
viewed as a random walk of the photon in the medium 
with step £ approximately. Thus the photon is multiply 
scattered about 6^ times before escaping. The Doppler 
shift is on average a random variable with zero average 
and about (fcw)^ variance. If successive scattering events 
are statistically independent, the photon frequency per- 
forms itself a random walk with step about kv, and 
the typical accumulated Doppler shift after N scatter- 
ing events is of the order of y/Nkv. As iV ~ 6^ = L"^ /l}), 
the key parameter is thus 



kf)V L 



(40) 



Knowing the expression of the eigenvalue, it is possi- 
ble to obtain analytically the behavior of the eigenvector 
close to resonance (i.e. for 5^0). By using a Taylor 
expansion of g^.s at second order in Eq. ([37]), we obtain 
the second derivative of the eigenvector at the resonant 
frequency given by 



yii/L,so 



1 



<5-+0 fcgu^ 



fo 

r 



3X2 



Using the previous results on the modal approach for 
the specific intensity for large systems at long times, we 
can easily deduce the expressions of the energy density 
and the outgoing fluxes close to resonance. Keeping only 
the lowest mode, the specific intensity writes 



/(z,/i,i,5) (^,(5)sin -z exp [soi] 



(41) 



Thus, the energy density is given by 

'TT 



C/(z,t,<5) = ^^(5)sin 

Ctr (O) 



L 



exp [sQt] 



Using the expression of the energy velocity, Eq. ([55]) . and 
the expression of the second derivative of the eigenvector 
near resonance, we obtain the spectrum of the energy 
density near the resonance in the form 



U(S) oc 1 

5^0 



1 



2k^v^ 



So 

r 



3L2 



This expression shows that it is impossible to have two 
peaks in the spectrum of the energy density for large 
systems at long times, because the second derivative re- 
mains negative whatever the size of the system. This 
is in complete agreement with the full Monte Carlo nu- 
merical simulations. To derive the same result for the 
outgoing fluxes, we have to integrate over the directional 
and the space variables the Fokker-Planck equation (i.e. 
Eq. ([34])) which gives 



+^ /Id 

fi[I{fi,L,S,t)~I{fi,0,S,t)]dfi^ (l+ \(3{S) 



Tdt 



-1 pL 



-I pL 



I (/i, z, 5, t) dzd/i + 

x/(/i',z,5,i)+7(<5) I' ifi',z,6,t) 
+1] (S) I" (/i', z, 5, t)] dzd/x'. 



(42) 



If ^ ^ 1, we are in the quadratic regime and the tem- 
poral decay rate r is given by the diffusion approxima- 
tion prediction (as for fixed atoms). This is nothing but 
the inequality ([3]) derived in the introduction using hand 
waving arguments. For ^ ^ 1, we are in the Doppler 
regime where the expression of the temporal decay rate 
is no longer quadratic in L. 



At long times, the entering flux vanishes. So the first 
term of Eq. (|42p is the total outgoing flux denoted by 
(j){t,S). This is the sum of the reflected {R{t,S)) and 
transmitted {T{t,S)) fluxes. Actually, at long times, a 
dynamical equilibrium is reached inside the system and 
the spectra of the reflected and transmitted fluxes are 
identical. Inserting the expression of the specific intensity 
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reduced to the lowest mode, Eq. pT]) . and using the dis- 
persion relation, Eq. p7p . the total outgoing flux writes 



APPENDIX A: SPATIO-TEMPORAL FOURIER 
TRANSFORMS OF RELEVANT QUANTITIES 



<^(t,(5) =2 . 



arctan 



^L{l + so/T)l3{6) 
L(l + so/r)/3(,5) 



(43) 



Using the same result as for the energy density, we ob- 
tain the second order approximation in frequency of the 
outgoing flux as 



fo 

r 



3X2 



Y2 



The sign of the second derivative of the flux may change 
and become positive for very large systems (i.e. L — > oo 
and So 0) which is at the root of the observation of 
two peaks on the reflected and transmitted fluxes. This 
is also in complete agreement with the numerical results. 
The physical picture is quite clear: far detuned photons 
are relatively rare in the medium, because they are less 
trapped and escape more rapidly: this is why they man- 
ifest themselves in the transmitted and reflected fluxes 
more strongly than in the energy density in the bulk. 



VI. CONCLUSION 

In this paper, we have derived a transport equation for 
the incoherent radiation propagating in a cold atomic gas 
at finite temperature. The derivation is based on first 
principles generalized to the case of moving scatterers. 
This equation, valid for the case of partial frequency re- 
distribution, is fully justified by microscopic arguments. 
It is solved numerically by an original and fully justi- 
fied Monte-Carlo scheme which gives reliable results on 
the outgoing fluxes and the energy density for all opti- 
cal thicknesses and incident beams. A modal approach 
is used to obtain information on the spectral and tempo- 
ral behaviors at long times. The main result is that the 
temporal behavior is a single exponential decay for all fre- 
quency components and the outgoing fluxes can exhibit 
two spectral peaks for very large systems. The velocity 
of energy propagation in the system is also affected by 
the motion of the atoms but is still reduced by typically 
four orders of magnitude compared to its velocity in vac- 
uum. This is an important result if we think of quantum 
memory applications for such systems. 
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Here is the list of all conventions used for the spatio- 
temporal Fourier transforms of the Green functions Gq, 
(G) , G (respectively in the vacuum, in the scattering sys- 
tem on average or not, all denoted by G in the following), 
of the mass operator Mm, of the intensity operator Km 
and of the scattering operator i^. 

G(r,i) = / G{k,uj)exp\ik-r- iujt]^^, 

M,n (r, t) ^ J (k, t^) exp \ik ■ r - iujt] ^ ^ , 

Km (ri,r3,r2,r4,ii,i3,i2,i4) = 

J Km (ki,k3,k2,k4,CJi,LJ3,CJ2,t^4) 

X exp [iki • ri - ik^ ■ ~ ik2 ■ V2 + ik^ ■ r^] 
X exp [—iujiti + ioJsts + iix'2^2 ~ iujit^] 
d^ki d^k2 d^ks d^k4 dwi du)2 dujs dw4 
87r3 87r3 Stt^ Stt^ '2n^~2n^' 

tm (ri,r2,ii,i2) = J t,n (ki,k2, wi, W2) 

X exp [iki ■ Vi — ik2 ■ r2 — iuJiti + iuj2t2\ 
d^ki d^k2 dwi dw2 
87r3 87r3 '2tt^' 

The signs in the exponentials are chosen such that the 
vertex operator Km describes correctly the correlation 
between the following two scattering processes: 



kz.ujz ~^ki,uji k4,tj4 k2,cJ2 



(Al) 



where k3 , k4 and ki , k2 are the incident and emergent 
wave- vectors respectively. W3,cij4 and ^^^1,^2 are the inci- 
dent and emergent frequencies respectively. 



APPENDIX B: MONTE CARLO SCHEME 



The temporal dependence of Eq. does not allow us 
to use a Monte-Carlo scheme to solve it numerically. Nev- 
ertheless, its spatio-temporal Fourier transform as writ- 
ten in Eq. ([27]) has the same structure as the standard 
RTE in the steady-state regime. This is the reason why 
we are using a Monte-Carlo type simulation. In prac- 
tice the iraplementation is the same as for the standard 
RTE [3^ |43| except that the probability densities are 
different and a Fourier transform is needed. Defining an 
effective extinction coefficient by 



^e.eff ((5, f^) 



fie {S - kou ■v,n)g (v) dv , 

Co 
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this equation writes 

/(q, u',n,S + ko (u'-u)-v) 



X 9 (v) 



dvdu'. (Bl) 



To go back to the real space for the space variable, we 
remark that 3? [fie,cS {S, ^)] is strictly positive for all 
and S which implies that 



1 



iq • U + He,cS 

and then 

1 



cxp [- (iq • u + ^e,cff) s] ds (B2) 



FT. 



iq • U + ^e.cff 



/•oo 

/ d {r - us) exp [~^e,cffs] ds 
Jo 

(B3) 

where FTq denotes the Fourier transform operator over 
q. Finally, the integral form of the transport equation 
writes 



To compute it, we use a Gauss-Hermite quadrature. 
Then, we perform the Monte-Carlo simulation The prob- 
ability densities used to compute the integrals are the 
following: 



the probability density to have an extinction pro- 
cess (in the sense of the coefficient /ie,Gff) at the po- 
sition s without having one from to s is ps (s) — 

He,eS (<5, ft) exp [-He,oS (S, fl) s] ; 



• the probability density for a photon to be scattered 
in the direction u coming from the direction u' is 
(u,u') = 1/(4^); 



Hr,u,n,5) 
1 



1 



^exff (S, fi) 



r+oc 

/ ^e,off (<5, exp [-^e,off ((5, f2) S] 
Js=0 



^ ^ j y ('^ ^ '^ou • V, VI) 

X / (r — su, u', fi, (5 + fco (u' — u) • v) dvdu'ds. 



(B4) 



This form fully justifies the use of a Monte-Carlo scheme. 
First, we have to compute the extinction coefficient for all 
frequencies fl and detunings 5. This coefficient is given 
by Eq. ^ (Voigt profile): 



/Xe,efl {6, ft) = 



4i7r/9 

2fcQ 



r/2 



-fcoU-v + f7/2 + ir/2 

r/2 



5 ~ko\i-w -n/2~iT/2 



1 



■ exp 



V 

2^ 



ift 

dv . 

Co 



• the probability density for an atom to have the ve- 
locity V is given by .g (v). Thus, the new frequency 
is given hy 6 = 5' — (u' — u) • v. 
The main difficulty is that the probability density ps is 
complex. In that case, we have to deal with the modulus 
of Ps as follows 



Ps [s) f {s)As 



„_,oo n ^ \ps (Sj)l 



(B5) 



where / 



C. The elements s,- are distributed in 



respect of the probability density \ps\. Note that the 
error bars on / (r, u, 17, 6) in the Monte-Carlo simulation 
have to be sufficiently small for the Fourier transform to 
give reliable results on / (r, u, t, 5). 
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